# Introduction ------
# Climate, Conflict, & Context: Re-Evaluating Americans' Support for Refugees
# Replication Materials
# Last Updated June 17, 2024
# Last Updated by Evan Stewart

# Loading Required Packages ------
library(tidyverse)
library(lubridate)
library(haven)
library(naniar)
library(infer)
library(cowplot)

# Load Data -------
pilot22<-read_dta("FILE PATH HERE")

# Data Cleaning and Recoding ------

clean<-pilot22 %>% 
  dplyr::select(starts_with("BBG6_"), Weights, partyid4, BBG5_13_3, ppethm, pppa1648, ppp20071) %>% 
  zap_labels() %>% 
  replace_with_na_all(condition = ~.x %in% c(-1, -2)) %>% 
  mutate(control = ifelse(is.na(BBG6_4), 0, 1),
         afghan = ifelse(is.na(BBG6_5), 0, 1),
         ukraine = ifelse(is.na(BBG6_6), 0, 1),
         climate = ifelse(is.na(BBG6_7), 0, 1),
         contact = recode(BBG5_13_3,
                          `1`=0,
                          `0`=1),
         white = ifelse(ppethm==1, "White","Non-White"),
         xtian = ifelse(pppa1648 %in% c(1:4,7,11), "Christian","Non-Christian"),
         borna = ifelse(ppp20071 ==1, "Born Again", "Not Born Again")) %>% 
  rowwise() %>% 
  mutate(DVadmit = sum(c_across(BBG6_4:BBG6_7), na.rm=T)) %>% 
  filter(DVadmit>0) %>% 
  mutate(DVadmit = 5-DVadmit)

clean$treat<-NA
clean$treat[clean$control==1]<-"Control"
clean$treat[clean$afghan==1]<-"Treatment:\nAfghanistan"
clean$treat[clean$ukraine==1]<-"Treatment:\nUkraine"
clean$treat[clean$climate==1]<-"Treatment:\nClimate"


# Table 2: Main Hypothesis Tests & Subsamples -------

# Main Effects
# NOTE: Main effect fails levene test for homogeneity of variance
# car::leveneTest(DVadmit~treat, data=clean)
# So, we relax the homogeneity of variance assumption using the oneway.test

clean %>% 
  select(DVadmit, treat) %>% 
  group_by(treat) %>% 
  summarize(out = mean(DVadmit, na.rm=T))

oneway.test(DVadmit~treat, data=clean)
chisq.test(table(clean$treat, clean$DVadmit))

## Pairwise t-test
pairwise.t.test(clean$DVadmit, clean$treat,
                p.adjust.method = "BH", pool.sd = FALSE)

# Republicans
clean %>% 
  filter(partyid4==1) %>% 
  select(DVadmit, treat) %>% 
  group_by(treat) %>% 
  summarize(out = mean(DVadmit, na.rm=T))

oneway.test(DVadmit~treat, data=filter(clean, partyid4==1))

clean %>% 
  filter(partyid4==1) %>% 
  mutate(DVadmit = as_factor(DVadmit)) %>% 
  chisq_test(DVadmit ~ treat)

# Democrats
clean %>% 
  filter(partyid4==2) %>% 
  select(DVadmit, treat) %>% 
  group_by(treat) %>% 
  summarize(out = mean(DVadmit, na.rm=T))

oneway.test(DVadmit~treat, data=filter(clean, partyid4==2))

clean %>% 
  filter(partyid4==2) %>% 
  mutate(DVadmit = as_factor(DVadmit)) %>% 
  chisq_test(DVadmit ~ treat)

# White
clean %>% 
  filter(white=="White") %>% 
  select(DVadmit, treat) %>% 
  group_by(treat) %>% 
  summarize(out = mean(DVadmit, na.rm=T))

oneway.test(DVadmit~treat, data=filter(clean, white=="White"))

clean %>% 
  filter(white=="White") %>% 
  mutate(DVadmit = as_factor(DVadmit)) %>% 
  chisq_test(DVadmit ~ treat)

# Nonwhite
clean %>% 
  filter(white=="Non-White") %>% 
  select(DVadmit, treat) %>% 
  group_by(treat) %>% 
  summarize(out = mean(DVadmit, na.rm=T))

oneway.test(DVadmit~treat, data=filter(clean, white=="Non-White"))

clean %>% 
  filter(white=="Non-White") %>% 
  mutate(DVadmit = as_factor(DVadmit)) %>% 
  chisq_test(DVadmit ~ treat)

# Christian
clean %>% 
  filter(xtian=="Christian") %>% 
  select(DVadmit, treat) %>% 
  group_by(treat) %>% 
  summarize(out = mean(DVadmit, na.rm=T))

oneway.test(DVadmit~treat, data=filter(clean, xtian=="Christian"))

clean %>% 
  filter(xtian=="Christian") %>% 
  mutate(DVadmit = as_factor(DVadmit)) %>% 
  chisq_test(DVadmit ~ treat)

# Non-Christian
clean %>% 
  filter(xtian=="Non-Christian") %>% 
  select(DVadmit, treat) %>% 
  group_by(treat) %>% 
  summarize(out = mean(DVadmit, na.rm=T))

oneway.test(DVadmit~treat, data=filter(clean, xtian=="Non-Christian"))

clean %>% 
  filter(xtian=="Non-Christian") %>% 
  mutate(DVadmit = as_factor(DVadmit)) %>% 
  chisq_test(DVadmit ~ treat)


# Figure 1: Average Support for Refugee Admission, Full Sample -----
# Make Plot
fig1<-clean %>% 
  group_by(treat) %>% 
  summarise(avg = mean(DVadmit),
            sd = sd(DVadmit),
            count = n()) %>% 
  mutate(se = sd/sqrt(count),
         lower = avg - (1.96*se),
         upper = avg + (1.96*se))%>% 
  ggplot(aes(x=treat, y=avg))+
  geom_point()+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=0.25)+
  geom_hline(yintercept=2.48, linetype=2)+
  scale_y_continuous(limits=c(1,4))+
  labs(x=NULL, y="Average Support for Admission")+
  theme_bw()
# Save as .TIFF
ggsave(filename = "~/Desktop/fig1.tiff", plot = fig1, dpi=300, width= 5, height = 3, units = "in")

# Figure 2: Average Support for Refugee Admission, Key Subgroups -----
p2a<-clean %>% 
  group_by(treat, partyid4) %>% 
  summarise(avg = mean(DVadmit),
            sd = sd(DVadmit),
            count = n()) %>% 
  mutate(se = sd/sqrt(count),
         lower = avg - (1.96*se),
         upper = avg + (1.96*se))%>% 
  filter(partyid4 %in% c(1,2)) %>% 
  mutate(party = ifelse(partyid4==1, "Republican","Democrat")) %>% 
  ggplot(aes(x=treat, y=avg))+
  geom_point()+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=0.25)+
  geom_hline(yintercept=2.48, linetype=2)+
  scale_y_continuous(limits=c(1,4))+
  labs(x=NULL, y=NULL)+
  theme_bw()+
  facet_wrap(vars(party))

p2b<-clean %>% 
  group_by(treat, white) %>% 
  summarise(avg = mean(DVadmit),
            sd = sd(DVadmit),
            count = n()) %>% 
  mutate(se = sd/sqrt(count),
         lower = avg - (1.96*se),
         upper = avg + (1.96*se))%>% 
  na.omit() %>% 
  #mutate(contact = ifelse(contact==1, "Close Contact","Non-Close Contact")) %>% 
  ggplot(aes(x=treat, y=avg))+
  geom_point()+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=0.25)+
  geom_hline(yintercept=2.48, linetype=2)+
  scale_y_continuous(limits=c(1,4))+
  labs(x=NULL, y=NULL)+
  theme_bw()+
  facet_wrap(vars(white))

p2c<-clean %>% 
  group_by(treat, xtian) %>% 
  summarise(avg = mean(DVadmit),
            sd = sd(DVadmit),
            count = n()) %>% 
  mutate(se = sd/sqrt(count),
         lower = avg - (1.96*se),
         upper = avg + (1.96*se))%>% 
  na.omit() %>% 
  #mutate(contact = ifelse(contact==1, "Close Contact","Non-Close Contact")) %>% 
  ggplot(aes(x=treat, y=avg))+
  geom_point()+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=0.25)+
  geom_hline(yintercept=2.48, linetype=2)+
  scale_y_continuous(limits=c(1,4))+
  labs(x=NULL, y=NULL)+
  theme_bw()+
  facet_wrap(vars(xtian))

# Make Plot
fig2<-plot_grid(p2a, p2b, p2c, ncol=1)

# Save as .TIFF
ggsave(filename = "~/Desktop/fig2.tiff", plot = fig2, dpi=300, width= 6, height = 7.5, units = "in")

# Figure 3: Google Trends ------

# Read in Data
g_refugee<-tibble::tribble(
  ~week, ~Afghanistan, ~Ukraine,
  "7/11/21",            1,      0.5,
  "7/18/21",          0.5,      0.5,
  "7/25/21",          0.5,      0.5,
  "8/1/21",          0.5,      0.5,
  "8/8/21",            3,      0.5,
  "8/15/21",           27,      0.5,
  "8/22/21",           15,      0.5,
  "8/29/21",           10,      0.5,
  "9/5/21",            3,      0.5,
  "9/12/21",            2,      0.5,
  "9/19/21",            1,      0.5,
  "9/26/21",            1,      0.5,
  "10/3/21",            1,      0.5,
  "10/10/21",            1,      0.5,
  "10/17/21",            1,      0.5,
  "10/24/21",            1,      0.5,
  "10/31/21",            1,      0.5,
  "11/7/21",            1,      0.5,
  "11/14/21",          0.5,      0.5,
  "11/21/21",          0.5,      0.5,
  "11/28/21",          0.5,        1,
  "12/5/21",          0.5,        2,
  "12/12/21",          0.5,        1,
  "12/19/21",          0.5,        1,
  "12/26/21",          0.5,        1,
  "1/2/22",          0.5,        1,
  "1/9/22",          0.5,        1,
  "1/16/22",          0.5,        3,
  "1/23/22",          0.5,        9,
  "1/30/22",          0.5,        4,
  "2/6/22",          0.5,        7,
  "2/13/22",          0.5,       15,
  "2/20/22",            1,      100,
  "2/27/22",            1,       89,
  "3/6/22",          0.5,       42,
  "3/13/22",          0.5,       29,
  "3/20/22",          0.5,       20,
  "3/27/22",          0.5,       13,
  "4/3/22",          0.5,       12,
  "4/10/22",          0.5,       11,
  "4/17/22",          0.5,       10,
  "4/24/22",          0.5,        8,
  "5/1/22",          0.5,        8,
  "5/8/22",          0.5,        8,
  "5/15/22",          0.5,        7,
  "5/22/22",          0.5,        6,
  "5/29/22",          0.5,        6,
  "6/5/22",          0.5,        5,
  "6/12/22",          0.5,        5,
  "6/19/22",            1,        4,
  "6/26/22",          0.5,        4,
  "7/3/22",          0.5,        4
) %>% 
  pivot_longer(cols=-week, names_to = "country", values_to = "search") %>% 
  mutate(week=mdy(week))

# Make Plot
fig3<-g_refugee %>% 
  ggplot(aes(x=week, y=search, linetype=country, group=country))+
  #geom_point()+
  geom_line()+
  scale_linetype_manual(name="Country",values=c(1,4))+
  #scale_color_manual(name="Country",values=c("gray50","black"))+
  labs(x=NULL,
       y="Normed Google Search Interest")+
  theme_bw()

# Make TIFF
ggsave(filename = "~/Desktop/fig3.tiff", plot = fig3, dpi=300, width= 6, height = 3, units = "in")
